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Abstract. Diffeomorphic matching (only one of several names for this technique) 
is a technique for non-rigid registration of curves and surfaces in which the curve or 
surface is embedded in the flow of a time-series of vector fields. One seeks the flow 
between two topologically-equivalent curves or surfaces which minimises some metric 
defined on the vector fields, i.e. the flow closest to the identity in some sense. 

In this paper, we describe a new particle-mesh discretisation for the evolution of 
the geodesic flow and the embedded shape. Particle-mesh algorithms are very natural 
for this problem because Lagrangian particles (particles moving with the flow) can 
represent the movement of the shape whereas the vector field is Eulerian and hence 
best represented on a static mesh. We explain the derivation of the method, and 
prove conservation properties: the discrete method has a set of conserved momenta 
corresponding to the particle-relabelling symmetry which converge to conserved 
quantities in the continuous problem. We also introduce a new discretisation for the 
geometric current matching condition of (Vaillant and Glaunes, 2005). We illustrate 
the method and the derived properties with numerical examples. 
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For Darryl Holm on the occasion of his 60th birthday 
1. Introduction 

Diffeomorphic matching is a numerical framework for quantifying the differences 
between geometric information (such as curves, surfaces, images or vector fields) using 
deformations from one geometric object to another. The geometric object is embedded 
in the flow generated by a time-series of vector fields with a chosen norm (such as 
the i7 1 -norm) defining the distance along the path between geometric objects. The 
computational task is to calculate the velocity fields which minimise this distance 
such that one geometric object is mapped to another. This framework was originally 
introduced in a series of papers including [GM981 ICYOll IMY014 IMTY03j being extended 
to match distributions (which can model curves and surfaces) in [GTY04j . and to 
match geometric currents (also for modelling curves and surfaces) in |VG05| . Various 
numerical approaches have been proposed for solving the optimisation problem, either 
by optimising the functional directly (with an extra term to penalise flows which do not 
map close to the target shape) as described in [BMTY05] , or by solving the equations of 
motion and shooting for a match between shapes by adjusting the initial conditions as in 
|TMT02[ IMM06at IMMS06j . The main challenge remains to find a numerical approach 
which is accurate and efficient, since the problem of computing the shortest path is a 
high-dimensional optimisation problem. 

In this paper we introduce a new numerical discretisation for the diffeomorphic 
matching problem in the context of matching (although it can also be used for matching 
surfaces, images and vector fields). This method uses a similar approach to the 
Hamiltonian Particle- Mesh (HPM) method |FGR02j . with the difference being that 
HPM uses the particle-mesh discretisation to interpolate density from the particles to the 
mesh, whereas in this application the particle-mesh discretisation is used to interpolate 
momentum (which takes a central role in the diffeomorphic matching framework). 

In section [2] we give a review of the diffeomorphic matching approach applied to 
curves in the plane, and establish the notation which will be used in the other sections. 
We also discuss the role of momentum, the conditions used to establish whether the 
curves have been matched, and the implications of the particle-relabelling symmetry 
satisfied by the equations of motion, as well as the connection with EPDiff and the 
Camassa-Holm equation. In section [3] we introduce the particle-mesh discretisation and 
discuss the discrete symmetries and conservation laws, as well as a discretisation of 
the current matching condition and a description of solution methods. In section [5] we 
illustrate the properties of the numerical method applied to computing the shortest path 
between two test shapes. Finally, in section [6] we give a summary and outlook. 
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2. Diffeomorphic matching of embedded curves 

In this section we describe the problem of matching one embedded curve onto another 
using diffeomorphisms. For simplicity we shall focus on simple closed curves in the plane 
although the approach is easily generalised to other structures. 

2.1. Parameterisations of curves 

We take two simple curves embedded in M 2 , Ca and Cb, with parameterisations 

Qa(s), Qb{s), sG[0,2tt). 

Note that the curves can equally well be represented by reparameterisations of the curves 
i.e. by 

Qa(v(s)), Qb(v(s)), sG[0,2tt), 

where t](s) is a diffeomorphism of the circle. Our aim is to find a matching process 
which is independent of such reparameterisations. 

2.2. Curves embedded in a flow 

If we take a vector field which defines a fluid flow u(x,t), and embed a curve in that 
flow, then the curve satisfies 

^Q(t;s)=u(Q(t;s),t), (1) 

i.e. each point of the curve moves with the vector evaluated at that point. 

The aim of the calculation is to search amongst time-series of vector fields u(x,t), 
t G [0, 1] such that ([I]) is satisfied, with the boundary conditions 

Q(0;s) = Q A (s), Q(l;s) = Q B (v(s)) 1 

for some (unspecified) reparameterisation t]. If these conditions are satisfied then we 
say that Qa is matched onto Qb by the vector field time series u(x,t). 

2.3. Optimisation problem 

We choose a norm for vector fields, such as the H™ norm defined by 

|| u ||| n = f u{x) ■ (1 - a 2 V 2 ) n u{x)d\o\{x) 
Jr 2 

Given a time series u(x, t), we can measure the total amount of deformation in the flow 
generated by u by the integral 

S(u) = jf^llufdi. (2) 

We wish to find the flow that maps Ca to Cb which is "nearest to the identity" with 
respect to the choice of norm. This leads to the following optimal control problem 
[GTY06] : 
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Definition 1 (Optimisation problem for curve matching) Let Qa, Qb '■ S 1 — > 

M. 2 be parameterisations of two curves Ca, Cb in the plane, and let || • || 2 be a norm for 
vector fields in the plane. Then the distance between curve Ca and Cb is defined to be 
the minimum over all vector fields u of the functional 

-\\u\\ 2 dt 



o 



2 

subject to the following constraints: 



• Dynamic constraint: ^Q(t] s) = u(Q(t; s),t), V£ G [0, 1]. 

• Matching conditions: Q(0; s) = Qa{s), Q(l; s) = Qb(t)( s ))> where r\ is some 
reparameterisation of the circle. 

2.4- Momentum 

We can enforce the dynamic constraint by introducing Lagrange multipliers P(t; s), so 
that the optimisation problem becomes 

8 [ ~\W\ 2 + [ P ■ (Q ~ u(Q))dsdt = 0, (3) 
Jo 2 J s=0 

subject to the matching conditions give above. The Euler-Lagrange equations then give 

P5(x - Q)ds, (4) 

ou J s=0 

Q =u(Q), (5) 
P = -P.Vu(Q), (6) 
where / = ||w|| 2 /2. For example, if we choose the if^-norm then 

= (1 - a 2 W 2 ) n u. 

ou 

We see that optimal velocity fields take the form 

/ PG(x - Q)ds, (7) 

where G is the Green's function associated with the chosen norm, e.g. for the if^-norm 
it is the Green's function for the operator (1 — a 2 V 2 ) n . 

Equation (j4j) allows us to write u as a function of P and Q, leading us to notice 
that equations ([MED are canonically Hamiltonian, with Hamiltonian function given by 

H = l(u(P,Q)) = ±\\u(P,Q)\\ 2 , 

i.e., half the square of the norm of the velocity field written as a function of P and 
Q. It is for this reason that we refer to the Lagrange multiplier P as the momentum 
associated with the curve. This Hamiltonian structure arises from the fact that the 
equations have been derived from a variational principle defined by the extreme points 
of the functional. Since the Hamiltonian for this system is time-independent, it is 
conserved along the trajectory, and hence the norm of the velocity is also conserved. 



The variational particle-mesh method 



5 



2.5. Matching condition 

There are a number of different possible ways to pose the matching condition 
mathematically. One matching condition which is invariant under reparameterisations 
is based on defining a singular vector field 

v Q (x) = r ^-5{x - Q)ds. (8) 

Js=0 OS 

We define a functional for these singular vector fields 
/ [ V Q] = J V Q -K* v Q dVo\(x) 

•£^ if(Q(s) - Q(s ' ))dsds ' (9) 

where K is some smooth kernel function. When Q matches Qb then f[v® — v® B ] 
vanishes. This is called the current matching condition |VG05] . 

This is a weaker condition that setting the value of Q(s) for each s at time t — 1, and 
hence we need to add another boundary condition to get a unique solution to equations 
(j4]l6]). In |MTY03j it was shown that the solution which minimises the functional is the 
one with P initially normal to the curve i.e. 

P . ?R. = 0, at t = 0, 

as 

and this extra boundary condition means that equations (jHl6]) have a unique solution. 

2.6. Relabelling 

The optimisation problem in definition [1] is invariant under symmetries given by 
reparameterisations of the curve 77. Noether's theorem tells us that the equations of 
motion (j4}|6]) have conserved quantities which are generated by these symmetries. 

To compute the conserved quantities we compute the infinitesimal generators of 
these symmetries 

5Q{s-t) = ^-i{s) 

where £ (s) is a vector field on the circle. The cotangent lift of this infinitesimal symmetry 
is 

5 (Q(s; t), P(s; t)) = ■ £(s), -P(s; t) ■ ^ • fa] 

which is a Hamiltonian flow generated by the Hamiltonian functional 

By Noether's theorem, is conserved along solutions of (jlj)-®, and since £(s) is 
arbitrary, we see that 

„ N dQ 
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is conserved along solution trajectories, for each s. In particular, this means that if the 
momentum is initially normal to the curve i. e. 

P( S ;0)-^(s;0) = 0, 

then the momentum is normal to the curve for all values of t along the solution. 
Therefore, optimal velocity fields take the form of equation (ED) with the added constraint 
that the momentum P is normal to the curve. 

2.7. EPDiff 

As described in |CH[ ICHH07] , because the dynamic constraint is given by a Lie algebra 
action of velocity fields on embedded curves, it is possible to eliminate P and Q by 
taking the time derivative of equation (jlj) and making use of equations ([ME])- This leads 
to a PDE defined on the whole of M 2 given by 

SI 

m, + V ■ (urn) + (Vu) m = 0, m = — . 

ou 

This is the Euler-Poincare equation on the diffeomorphism group, abbreviated as EPDiff 
[HM04t IHMR98j . which is the equation for geodesic flow on the diffeomorphism group. 
Although we do not explicitly solve EPDiff during the computation of the optimal flow, 
it is useful to understand the computed solutions as singular solutions of EPDiff. In one 
dimension and with the if^-norm, EPDiff becomes the Camassa-Holm equation |CH93] 
which is completely integrable with singular soliton solutions. 

3. Particle-mesh discretisation 

We take a discrete mechanics and optimal control |JMOB05] approach by applying the 
discretisation directly to the functional ([3]) and deriving the resulting equations. The 
discretisation used is a particle-mesh method with: 

• the vector fields discretised on a fixed, finite mesh, and 

• the curve discretised set of moving points. 

The principle benefits of this discretisation are that different norms can be defined on 
the mesh without the need to calculate Green's functions. With large numbers of points 
the method becomes very efficient as particle momenta can be interpolated to the mesh, 
then the norm operator can be inverted, then the velocity values can be interpolated 
back to the particle positions. 

3.1. Mesh discretisation of vector fields 

We take a fixed set of points on a mesh {x^^Li (f° r the numerical examples in this 
paper we used an equispaced square mesh) and give each mesh point x k a vector u k . 
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We interpolate from the set {wfe}^ of vectors to a general point x in the plane using 
a linear interpolation 

Tig 

u(x) = ^2u k ifj k (x). (10) 

k=l 

For the examples in this paper we set ipk to be a tensor product of cubic B-spline 
functions centred on x = xj,. 

Remark 2 The set of vectors on the mesh generate a finite dimensional subspace of 
the infinite dimensional space of vector fields. However, the finite dimensional subspace 
is not closed under the Lie bracket for vector fields which means we will not be able to 
obtain a discrete EPDiff equation from the discrete equations of motion by eliminating 



Q and P as in section {2. 7| ). 



Once we have this discrete representation of the vector fields we can define a 
discretised norm using standard mesh methods (such as finite difference, finite volume 
etc.). For the the examples in this paper we took periodic boundary conditions, and used 
a spectral discretisation of the (1 — a 2 V 2 ) n operator using discrete Fourier transforms, 
and a simple Riemann sum for the integration (which gives spectral accuracy in this 
case). Since the Green's functions decay exponentially over the lengthscale a, the 
boundary conditions should not affect the solution as long as the curve is sufficiently 
far from the boundary. Other boundary conditions are possible if, for example, a finite- 
difference method is used to discretise the norm operator. 

3.2. Particle discretisation of curves 

We replace the parameterised curve Q(s) by a finite set of points in the plane {QpYaLi- 
The equation (0Q) gets replaced by the semi-discrete (continuous time/discrete space) 
equation 

Tig 

Qp = ^ u k4>k{Q(3)- (ii) 

k=l 

We also have to choose a discretisation of integration around the loop; again for the 
examples in this paper we use a Riemann sum. 

3. 3. Semi-discrete functional and equations of motion 

After the particle-mesh discretisation we obtain a semi-discrete functional by integrating 
the discrete vector field norm from t = to t — 1 and introducing Lagrange multipliers 
to enforce the constraint (TTTT) to get 

I \\H\ 2 9 + Y. p v (^Qe-^kMQ^dt. (12) 
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The Euler-Lagrange equations for the optimal values of this functional are 

\\\< = E^(^). ( 13 ) 



du k 2 

p 

n g 

Q/3 = ^U k ^ k (Qp), (14) 

k=l 

Pfs = -Pp-Y, u ^MQp)- (15) 

k 

Once again, equation ( fl3i) allows us to write u as a function of P and Q, leading us to 
notice that equations fll4H15p are canonically Hamiltonian, with Hamiltonian function 
given by ||u(P,Q)||J/2. 

3.4- Time discretisation 

Equation ffTTj) can be discretised using any one-step method, such as a Runge-Kutta 
method. For simplicity, in this paper we use the forward Euler discretisation, but the 
derivation of the equations is very similar for other methods. The equation becomes 

Q; +1 = Q n p + At^ n+1 MQ n )- 

k 

The discrete functional becomes 

E ( a 4ii«ii; + E p p +1 ■ (op 1 - Q£ - a * E < +1 MQ n P ) ) ) , (is) 

«=1 \ P \ k=l J J 

where At = 1/N, and the discrete Euler-Lagrange equations are 

^\\^ n+ x=J2 p r 1 MQ^ (it) 

=Q- + AtJ2< +1 MQp), (18) 

k 

p; +1 = p; - Atp n+1 ■ K +1 vMQp)- (is) 

k 

These equations provide a variational integrator |LMOW03] for the semi-discrete 
equations ( fT3lfT5l) . Hence, they give a symplectic integrator (see |LR05] for a survey of 
these methods) for the semi-discrete equations in Hamiltonian form. In this particular 
case we obtain the first-order symplectic Euler method; higher-order partitioned Runge- 
Kutta methods can be obtained by discretising ffTTj) using any Runge-Kutta method. 



3.5. Discrete symmetries 

In this section we discuss the properties of applying a variational integrator to equations 
(113j) - fll5j) . and their possible benefits for the problem of matching curves. 

The properties that make variational integrators the best choice for long 
integrations are: 
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• Modified Hamiltonian: if the Hamiltonian H(P, Q) is analytic then it is possible 
to use backward error analysis to find a modified Hamiltonian 



See |LR05j for a summary and references. 
• Conserved momenta: if the continuous-time Euler-Lagrange equations have 
conserved momenta associated with symmetries of the Lagrangian, then discrete 
Euler-Lagrange equations also conserve these momenta [LMOW03] . 

3.6. Modified Hamiltonian 

In the problem of matching curves, the Hamiltonian is half the squared norm for vector 
fields, which is conserved along trajectories as described in section 12.41 If a symplectic 
integrator is used then backward error analysis guarantees that the Hamiltonian will 
be approximately conserved for long times (although it should be noted that no theory 
exists for the piecewise-cubic functions used in the examples of this paper). It remains 
an open question as to whether the conservation of the Hamiltonian is important for 
these problems on finite-time intervals, since they are on short time intervals and so 
the value of the Hamiltonian will also be approximately preserved by variable step-size 
methods using error control. 

The main cost of using variational integrators is that they are implicit for this type 
of system where the Hamiltonian is a function of P and Q, but during the optimisation 
algorithm we are able to store values along the trajectory from previous integrations 
which can be used as initial guesses for solving the implicit equations using Newton 
iteration. The variational integrator allows large timesteps to be taken in that case 
(although of course this is the case for any implicit method, symplectic or otherwise). 

3. 7. Momentum conservation 

As described in section 12.61 the equations (EH) have a symmetry under 
reparameterisation of the curve which has an associated conserved momentum P ■ 
dQ I ds. This means that for the optimal solution, P remains normal to the curve along 
the trajectory. Since we have discretised the curve we have broken that symmetry, but 
as we shall show in this section, it is still possible to recover a sense in which P is normal 
to the curve along the solution. 

Equation ffTOl) defines a vector field over the whole space parameterised by the 
vectors {uk\1 9 =1 associated with the mesh points. This means that we can choose a 
parameterised curve Q{s) which passes through our discrete points {QpYpLx in sequence 
and follow its evolution along the flow by solving 



H = H + At p AH( 



P,Q;At) 



which is conserved over times 



t\ < Cq exp 



a /At 




(20) 
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More generally, we can follow how the flow generated by the vector fields evolves over 
the whole space 
d 

—g{x,t) = ^2u k ijj k (g(x,t)), 

k 

where g(x, t) is the flow map taking points from their position at time to their position 
at time t (x plays the role of a coordinate on "label space" here), and we can follow the 
Jacobian of this map 

d dg(x,t) , dg(x,t) 

k 

In particular, we can evaluate this equation at each of our discrete points Qp\ 

j p = Y,UkVMQf3)-J^ Jp = ^(Qt>(o),t), Mo) = id. (21) 

k 

A variation 5Qp(0) in the initial conditions Q/?(0) leads to a variation in the entire 
trajectory given by 

SQfi(t) = Jp{t)5Q p {$). (22) 

This variation generates a symmetry of the equations, as shown in the following lemma. 

Lemma 3 The infinitesimal transformation given by equations $%T\}-$EBj) together with 

5P[3 = 0, (3 = l,...,n p , 5u k = 0, k=l,...,n k , 
is a symmetry of equations [T3\) - [Th*\) . 

Proof. Since the equations have been derived from an action principle, we simply need 
to show that the infinitesimal transformation causes the action ffT2l) to vanish. If we 
apply the transformation to the integrand (the Lagrangian) then we obtain 



5L = 5 ( ~\\u\\ 2 g + P P ■ ( Qp ~ jt u ^MQf: 
\ p \ fc=i 



v > 



fc=i 



= Y J p p-[jp- 5 QM -J2^MQp) ■ JpSQp(o) 
p \ k=i 

(n g 
J2 u ^MQp) ■ Jp ■ SQp(o) -J2 u ^MQp) ■ Jp$Qp(o) 
k k=i 

= o, 

and hence the result. 
□ 

This symmetry of the equations has an associated conserved momenta following 
Noether's theorem, as described in the following proposition. 
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Proposition 4 For each /3 = 1, . . . ,n p , the quantity JjPp is conserved. 

Proof. Applying the symmetry to the action principle and substituting the equations of 
motion ([TBl) - (fToT) gives 



9 




p \ k=i 

^ k=l 




We obtain the result since 5Qp(0) is an arbitrary vector. 
□ 

As described in |LMOW03] . these conservation laws satisfied by the semi-discrete 
equations will be preserved by numerical methods provided that they are derived from 
a discrete variational principle i. e. following the framework described in section 13.41 

If we choose a continuous curve through our set of points that moves with the flow 
according to equation (1201) . then an approximation to dQ(sp)/ds ■ ds is given at time t 
by 

JpAQpiO), (3 = l,...,n p , 

where AQ j a(0) is an approximation to dQ(sp)/ds ■ ds at time t = 0. This leads to the 
following corollary: 

Corollary 5 If Pp is chosen to be normal to AQp(0) at time t = for all f3 i.e. P is 
initially normal to the curve, then P is normal to J / gAQ ( g(0) for all t along the flow, 
i.e. P stays normal to the curve. 

Proof. For each (3 = 1,..., n p , the component of P tangential to the shape is 
Pp(t) ■ Jp{t)AQp{0) = P p (0) ■ J^(0)AQ^(0) 
= P p (0) ■ IdAQ^(O) 
= P (O)AQ (O) = 0. 

□ 

If the equations of motion for Q and P are discretised using a variational integrator 
then this property is preserved, provided that the discrete equation for the evolution of 
J is defined as the gradient of the evolution equation for Q, since then it generates a 
symmetry of the discrete variational principle just as in the continuous time case. 

For example, the discrete equation for the gradient of the time-discrete flow for Q 
obtained from the symplectic Euler method is 

= r p + At Ys< +1 ^m j l ( 23 ) 
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Proposition 6 For each (3 = l,...,n p , the quantity JpPp is conserved when the 
equations are integrated using the symplectic Euler method, and J is obtained from 
equation \23\ 



Proof. 



N / / n g 

s a = 5 £ -At||«ii; + y, P T X ■ QT 1 - Q£ - A * E < +1 MQ} 

71=1 \ /3 \ k = l 

N 

= E ( E p ; +1 ■ ( tQT 1 - 8Q? - E < +1 v^(Q^) • ^ 

n.=l \ /3 \ fe=l 

^9 



N 



= E E p p +1 ■ ( J T l - J £) • - E < +1 vMQp) ■ JpSQj 

n=l /3 \ fe=l 

= E p « +1 • ( E < +1 vMQp) ■ n ■ SQl - J2 K +1 ^MQp) ■ JpSQ 

/3 \ k k=l 

= 0, 

Applying the symmetry to the discrete action principle [16] and substituting the equations 
of motion (Tl3|) - (fr5|) gives 

o = E ( p ; +1 J T l - p m) ■ 6 Q% 

We obtain the result since <5Q° is an arbitrary vector. □ 

We can verify the discrete conservation of PpJp directly since 

ppi = p; +1 (h + J2 u k +1 vMQp))Jp 

k 

pn+l jn+1 

- ^p Jp ■ 

Similar results can be obtained when higher-order variational integrators are used by 
discretising the Q equation using a Runge-Kutta method. 

As discussed in section 12.51 the norm-minimising solution has momentum normal 
to the shape along the whole trajectory. Good preservation of the tangential component 
of momentum is important because it means that it is only necessary to constrain the 
momentum to be normal to the shape in the initial conditions when a shooting algorithm 
is used (discussed in section @J. 

3.8. Discrete matching condition 

We can also apply a particle-mesh discretisation to the matching condition described in 
section 12.51 The particle-mesh representation of equation (JE} is 
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where AQ@ is an approximation to dQ(s)/ds ds at s = sp. For example, a simple finite 
difference approximation gives 

v ? = ^{Qp - Qf3-i)ipk(Qp)- 

This is the approximation which we will use in the computed numerical examples used 
in this paper. 

Once the singular vector field is evaluated on the mesh we can apply standard mesh 
discretisation methods to compute the functional ([9]): 

kl 

where K^i is a discretisation of the kernel operator K. For the computed numerical 
examples used in this paper, the kernel operator used was the inverse operator 
(1 — ct 2 V 2 ) -2 discretised using discrete Fourier transforms. 

This results in a numerical discretisation of the functional used for the current 
matching condition. After numerical discretisation the functional will not have a 
minimum at zero any more, and we must aim to minimise this functional rather than 
find the zero to achieve matching. 

3.9. Efficiency 

In |VG05] . a mesh-free method was applied to the current matching problem. The 
main cost in this type of method is in summing up Green's functions on each of the 
points on the curve to calculate their velocity, and the value of the matching condition. 
They employed a fast multipoles algorithm [GS91] which has a computational cost of 
size 0(N log N) (where N is the total number of particles), although the multiplicative 
constant in the scaling can be relatively big, requiring larger N before benefits are seen. 

For a particle-mesh method, the cost of evaluating the momentum on the grid is 
O(N), and the cost of the FFT is 0(M 2 log M) (where M is the number of rows of 
gridpoints i.e. the total number of grid points is M 2 ), although one could reduce this 
by discretising the operator used in the norm for velocity (e.g. the Helmholtz operator in 
the examples given here) and just performing a few iterations of a method such as Jacobi, 
resulting; in an 0(M 2 ) cost. In the case where one is matching dense information (images 
or measures, for example), then typically N = cM 2 (with c > 1) and the particle-mesh 
approach produces a method which is competitive with fast multipoles methods. For 
the case of matching curves, typically N = cM (with c > 1) and a good fast multipole 
implementation will be more efficient. However, any mesh-based operator inversion is 
easily parallelised using standard methods (as is the particle-mesh operation), and so 
this could be a strength of the particle-mesh approach in the curve case. There are also 
other advantages to using a mesh, for example it is very easy to modify the operator to 
investigate the behaviour with different kernels. 
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4. Solution methods 

As discussed in section 13.81 numerical discretisation of the matching condition means 
that we cannot require that the matching functional vanishes, and so we must take this 
into account in the solution approach. There are two main approaches to obtaining a 
numerical approximation to the optimal flow between embedded curves: 

(i) Inexact matching: In this approach, suggested in |MTY03| . we "softly" enforce 
the matching condition by adding a penalty term to the functional (j2J): 



One then attempts to minimise this functional directly by applying a descent 
algorithm (such as nonlinear conjugate gradients) by varying the time series of 
vector fields and computing the implicit gradient of the functional /. An alternative 
method is to introduce the dynamical constraint using the Lagrange multipliers 
P as in section 12.41 and to solve the resulting Euler-Lagrange equations (jMHl) 
(modified to accomodate the penalty functional) using Newton iteration. This 
approach becomes attractive when it is possible to find a good initial guess at the 
solution e.g. by specifying a path of curves between Ca and Cb and approximately 
solving equation (jSJ) to get an initial guess for P given the initial guess for Q. 
It is worth mentioning that the minimisation of the functional [24] becomes ill- 
conditioned as a — > and so it is not always possible to match one curve onto 
another with the desired accuracy (or to within the size of numerical errors, having 
discretised the functional). A modification of the functional called method of 
multipliers, described in |Ber82j . introduces Lagrange multiplier variables along 
with the penalty parameter a and allows a reduction of the error in matching 
without making a arbitrarily large. 

(ii) Minimisation by shooting: In this approach, advocated in [MM06b] . we seek 
initial conditions for P such that the functional [9] is minimised. The gradient of the 
functional with respect to the initial conditions for P is computed by solving the 
Euler-Lagrange equations (HKSD from t = to t = 1 with those initial conditions, 
computing the gradients of f[u QN ] with respect to Qp, and propagating those 
gradients back to P^ using the adjoint equations (see |Gun03] . for example). The 
problem can then be solved using a nonlinear gradient algorithm such as nonlinear 
conjugate gradients. As described in section 12.51 P must be constrained to be 
normal to the curve in order to obtain the optimal path. 
This is the approach which we used in computing the numerical examples. 

5. Numerical examples 

In this section we show a computed curve-matching calculation of two curves in the 
plane. The curve was parameterised using 420 points, with a square 128 x 128 mesh of 
size 27r x 2tt. The velocity norm used was the if^-norm with a = 0.4, discretised on 
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Figure 1. Test curves used for the example matching computation. 




Figure 2. Snapshots of the optimal path between two test curves computed using the 
particle-mesh discretisation. The computed curve is plotted with a continuous line, 
superimposed on the target curve plotted with a dashed line. Top row: the curve at 
t = [0, 0.2, 0.4]. Bottom row: the curve at t = [0.6, 0.8, 1]. 



the mesh using discrete Fourier transform, and the kernel used for the current matching 
was the Green's function of the (1 — a 2 V 2 ) -2 operator with a = 0.4. 
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Figure 3. Plots illustrating the optimal deformation map which maps between our 
two test curves. The numerical solution specifies a vector field which is defined 
everywhere, which can be used to transport other points which are not on the curve. 
This calculation was performed on a set of points on equispaccd grid lines to show 
how space is being deformed around the curve. Top: the grid lines and curve before 
deformation. Bottom: the grid lines and curve after deformation. Note that the flow 
vector fields are far from divergence- free, as can be seen by inspecting the areas of the 
squares on the deformed grid. Note also that the optimal flow only deforms space near 
to the curve. 
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Figure 4. Plots showing the evolution of the momentum P along the optimal flow. 
The plots are taken from the same snapshots as in figure [2l with vectors showing the 
direction and magnitude of P on the curve. Note that P remains normal to the curve 
throughout. 




Figure 5. Plots showing JpAQp around the curve which is an approximation to 
dQ/dsds as described in section [5751 This approximation remains normal to the curve 
throughout, with the increase in magnitudes showing regions of stretching in the flow. 
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The solution was obtained by the "minimisation by shooting" method as described 
in section HI The minimisation was performed using the Scientific Python [JOP + ] 
optimize . fmin_ncg routine which applies the Newton conjugate gradients algorithm 
with the hessian computed by applying finite differences to the gradient. The algorithm 
was run out until the functional was reduced to a value of 1.748 x 10~ 5 (with an initial 
value of 0.0108). 

A plot of the two curves used for the tests is given in figure [TJ These two curves 
are quite different and require large deformations to transform one curve into the other. 
The calculated path is illustrated in figure [2] with a few snapshots of the curve during 
the transformation at various times. The effect of the deforming flow on the surrounding 
space is illustrated in figure [3j 

To illustrate the results of section [331 a plot showing momentum vectors at various 
times is given in figure 0] which suggests that the momentum stays normal to the 
curve. This is confirmed by figure showing the evolution of J / gAQ / g(0) which is 
an approximation to dQ/dsds. The quantity JT^-Ps ■ JpAQp remains within round-off 
error of zero throughout, confirming the results of section [3751 

6. Summary and outlook 

In this paper we introduced a new particle-mesh discretisation for diffeomorphic 
matching of curves, which can also be used for matching images. In this method the 
vector fields used to transport the curves are represented on a fixed mesh, whilst the 
curves themselves are represented by a finite set of moving points. Since the discrete 
equations arise from discretising an action principle, they are variational integrators 
which have many favourable properties. As discussed in |MM06aj . whilst the benefits 
of variational integrators have been established for long integrations such as those for 
celestial mechanics or molecular dynamics, the benefits for short time optimal control 
problems such as the curve matching problem discussed in this paper are not so clear. 
There are a number of drawbacks and benefits which need to be investigated with 
further testing. In this paper we showed that the variational integrators that arise from 
the particle-mesh discretisation have a discrete form of the momentum conservation 
law which leads to the curve momentum remaining normal to the curve throughout 
the computed trajectory between shapes. We also noted that the integrators have a 
modified Hamiltonian which is conserved over long times (but not exponentially long 
since the compactly supported basis functions ipk{x) are not analytic) which can be 
interpreted as a modified metric for the discrete equations. We then showed illustrative 
examples obtained from the calculation of the optimal trajectory between two closed 
curves in the plane. 

The main computational challenge for diffeomorphic matching remains the design of 
efficient algorithms to obtain optimal trajectories for large datasets. Of the two solution 
methods described in section [H the inexact matching method applied to the particle- 
mesh discretisation may be best applied in parallel by using a Newton-Krylov method 
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to solve the Euler-Lagrange equations (including the penalty term) since it is possible to 
obtain a good initial guess for the optimal path by other methods. The minimisation-by- 
shooting approach might be best applied using a multilevel scheme where the problem is 
first solved with a small number of particles, with more particles being introduced once 
the lower dimensional problem has been solved to sufficient accuracy. These approaches 
will need to be developed in order for the discretisation to be applied to practical 
engineering applications, as well as convergence studies and error estimates. 

In future work we shall investigate the convergence properties of this method in the 
limit as the number of points in the discretisation of the curve goes to infinity (together 
with the number of mesh points), and apply the method to investigate practical datasets. 

Although we do not compute examples here, this method could also be used 
for matching surfaces in three dimensions using a current-matching approach. The 
main modifications are that equation ( flOl) needs to be evaluated in three dimensions 
using a tensor product of three B-splines (one for each Cartesian component), and the 
parameterisation of the surface should become a triangulation with a singular current 
being interpolated from the surface to a three-dimensional mesh using the same basis 
functions. We will develop and investigate such a method in future work. We will also 
investigate efficiency of the particle-mesh versus the mesh-free fast multipoles approach 
of |VG05] in numerical tests on real data. 

Further developments will be to apply the particle mesh to the related problem 
of metamorphosis [TY04] , and to problems where quantities such as vector and tensor 
fields associated with images also need to be matched together. 
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